EFFICIENT CALCULATION OF DETERMINANTS OF 
SYMBOLIC MATRICES WITH MANY VARIABLES 
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^_j Abstract. EfRcicnt matrix determinant calculations have been 

^-s studied since the 19th century. Computers expand the range of de- 

^sj terminants that are practically calculable to include matrices with 

» . symbolic entries. However, the fastest determinant algorithms for 

Oh numerical matrices are often not the fastest for symbolic matrices 

■^r with many variables. We compare the performance of two algo- 

rithms, fraction-free Gaussian elimination and minor expansion, 
on symbolic matrices with many variables. We show that, under a 
simplified theoretical model, minor expansion is faster in most sit- 
uations. We then propose optimizations for minor expansion and 
demonstrate their effectiveness with empirical data. 
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,__! 1. Introduction 

,__! Determinants of square matrices are essential in a myriad of fields, 

On both theoretical and applied. Computers can be used to quickly calcu- 

late determinants of matrices with numeric and, in more recent decades, 
symbolic entries. Efficient algorithms exist for matrices with a small 
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(^ number of variables of high degree. We turn our attention here to 

^^ the opposite case: matrices of polynomials with many variables of low 

degree. 

Bareiss introduced an algorithm, fraction-free Gaussian elimination, 
^ that requires 0{n^) polynomial operations and avoids fractions [Ij. 

j^ However, the cost of those operations becomes prohibitively large. An- 

other algorithm, minor expansion, takes 0(2"n) polynomial operations 
and has been shown to be an attractive alternative ^. Our first re- 
sult, presented in Section [3| is a quantitative comparison of the cost 
of these two algorithms on a dense matrix of linear polynomials. We 
find that minor expansion is favorable unless the size of the matrix is 
greater than some function / of the number of variables s, which by 
our evidence grows faster than linearly. To the extent that our com- 
puter hardware can handle, the theoretical analysis appears reasonable 
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when compared to experiments, though the exact correlation depends 
heavily on implementation details of each algorithm. 

Griss introduced row-ordering as a technique to reduce further the 
cost of minor expansion [3l H]. Our second result, presented in Sec- 
tion |4| suggests the scenarios when row sorting is most helpful, namely, 
in our experiment, when approximately half of the entries of a matrix 
are 0. We try several possible sorting strategies, all of which perform 
similarly, so further analysis is needed to determine if any of them have 
an advantage over the others. 

2. Preliminaries 

We use the following notation: 

A = a.n n X n matrix with entries in Z[x]. 
Qij = entry of A in row i and column j; starting index is 1. 
det{A) = determinant of A. 
[a] = {neZ:l<n<a}. 
suh{A, I, J) = submatrix of A using rows in / and columns in J, 
where I, J ^ [n]. 
T{p) = number of terms in polynomial p. 
Let i G [n] . The determinant of A can be defined recursively as 

J ail if n = 1 

1 Sfci(~l)"^"''^«i det(y4jj) otherwise, 

where Aij = sub(y4, [n] \ {i}, [n] \ {j}). Note that any choice of i 
yields the same result, as does summing over i with fixed j. Naively 
calculating det(y4) requires computing n! products of n factors each, 
which would take 0{{n + 1)!) polynomial multiplications. There are 
several algorithms that improve drastically on this performance, two of 
which we study here. 

2.1. Minor expansion. Calculating the determinant of A requires 
calculating the determinant of Aij = sub(y4, [n] \{i}, [n] \ {j}) for fixed 
i and all j G [n] . Similarly, the calculation of each of these determinants 
requires calculating the determinants of sub(A, [n] \ {i, k}, [n] \ {j, /}) 
for fixed distinct i and k and all distinct j,l G [n]. Notice that, for 
given j and /, the determinants of both Aij and An require calculating 
the determinant of sub(A, [n] \ {i, k}, [n] \ {j, /}). Taking advantage of 
this and similar redundancies allows us to reduce the number of deter- 
minants we have to calculate. Increasing i from 2 to n, we calculate 
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determinants of every i x i submatrix contained in the first i rows as 
a linear combination of the (i — 1) x (i — 1) determinants calculated 
in the step before. The determinants of these submatrices are called 
minors. The minor expansion algorithm follows. 

input A 

M0:=1 

% Calculate determinants of submatrices of increasing 

size. 

for i := 2, . . . ,n do 

for J C [n] : \J\ = i do 

% Calculate det(sub(yl, [i], J)). 

% J is sorted with k^^ element jk. 

endfor 
endfor 

% A has all n columns, so det{A) = M[„]. 
return M[„] 

Minor expansion requires X]r=2(0(") — 0(2"ra) polynomial multipli- 
cations, which, we will see, is more than fraction-free Gaussian elimi- 
nation requires. 

Not every entry of A is involved in the same number of multiplica- 
tions. Entries in the i^^ row are directly present in (^) multiplications 
except for those in the first row, which are in (2) multiplications. En- 
tries in earlier rows, through their presence in minors, are indirectly 
involved in many more multiplications. 



2.2. Fraction-free Gaussian elimination. It is clear from (2.1) that 
the determinant of a triangular matrix is the product of its diagonal 
entries. Ordinary Gaussian elimination, described by the iteration be- 
low, calculates a sequence of matrices A^^-*, starting with A^^' = A, such 
that A^"-^ is upper triangular and det{A^''^) = det(A) for all /c G [n]. 

(fc) (fc) 
4^'^ = oi^-^^ V(z,,)e{A: + l,...,n}x{fc,...,n}. 

Unfortunately, this method requires fractions, which, especially when 
working with polynomials, are costly to reduce with GCD calculations 
and costly to let grow without reduction. 

Fraction-free Gaussian elimination calculates a different sequence of 
matrices A^''\ again with A^^'^ = A, such that ann = det(y4). We define 
agj = 1. The one-step fraction-free Caussian elimination algorithm 
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calculates the sequence A^'^l with the iteration 

[k] [k] _ [k] [k] 

a\j J = 'jj—^^ '- y{i,j) e{k + l,...,n}x{k,...,n}. 

'^k-l,k-l 

It can be shown with Sylvester's Identity that the division has no re- 
mainder and that the divisor is the largest possible that guarantees 
this [Tj. The algorithm requires Y17=i 3(^ — 0^ — 0{n^) polynomial 
multiplications and divisions. 

From here on, "Gaussian elimination" will refer to one-step fraction- 
free Gaussian elimination unless otherwise specified. We note that 
there is a "two-step" variety of fraction-free Gaussian elimination, which 
calculates A''^' for only odd k using a more complicated formula, and 
that larger step sizes are possible. 

3. Algorithm comparison 

We mentioned previously that minor expansion and Gaussian elim- 
ination take 0(2"n) and 0(?2^) polynomial operations, respectively. 
However, the time taken by an individual polynomial operation can 
vary drastically depending on the polynomial, and minor expansion 
very often outperforms Gaussian elimination for a variety of reasons. 

3.1. Theoretical example. We first consider a model in which we 
examine a particular class of matrices and measure the cost of each 
algorithm in terms of integer operations. We make the following sim- 
plifying assumptions: 

• Addition and subtraction of polynomials have zero cost. 

• Multiplication and division of polynomials p and q take T{p)T{q) 
integer operations. 

• All integer operations have the same cost. 

Choose some s G N. Suppose that every entry of A is a linear 
polynomial in s variables where each term has exactly one variable as 
a factor. That is, aij = J2k=i^ijk^k for all i,j G [n], where all Xk are 
variables and all Cijk G Z. 

The chance of having zeros as coefficients of polynomials in the fol- 
lowing calculations is negligible, so we ignore the possibility for sim- 
plicity. 

Call a polynomial i-homogenous iff it has total degree i in each of 
its terms. For example, the entries of A are 1-homogenous. An i- 
homogenous polynomial has at most (*^1^ ) terms. Note that the prod- 
uct of an i-homogenous polynomial and a j-homogenous polynomial is 
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Figure 3.1. The ratio of the cost of minor expansion 
to the cost of Gaussian ehmination. Note that the origin 
is on the farthest vertical edge of the surrounding box. 



(i+j)-honiogenous and that the sum of two i-homogenous polynomials 
is also z-homogenous. 

We turn our attention first to minor expansion. When calculating 
the ixi minors of A, we perform i (") multiplications of a 1-homogenous 
polynomial by an {i — l)-homogenous polynomial. The maximum cost 
of minor expansion in integer operations, denoted Cm, is 
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The cost of Gaussian elimination can be computed similarly. When 
dealing with a particular entry in the i^^ elimination step, we do two 
multiplications of two z-homogenous polynomials and one division of a 
2i-homogenous polynomial by an {i — l)-homogenous polynomial. This 
happens to the lower {n — i) x {n—i) block of the matrix. The maximum 
cost of Gaussian elimination in integer operations, denoted Cg, is 
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Figure 3.2. Predicted crossover points between mi- 
nor expansion and Gaussian elimination versus observed 
crossover points. Matrices are n x tt, with 1-homogenous 
linear polynomials in s variables for entries. 
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3.2. Experimental boundary points. Figure 3.1 

with respect to n, the dimension of A, and s, the 

ables. Unless n grows much faster than s, minor expansion has less 

costly asymptotic behavior. In the range we calculated, n must grow 

faster than linear with respect to s for Gaussian elimination to be the 

better choice asymptotically. 

To test the predicted boundary points — points where the two algo- 
rithms have the same cost — we ran the following experiment. Starting 
from n = \ and s = 1, we randomly generate A and time the two 
determinant algorithms applied to it. When minor expansion is faster, 



we increment n; otherwise, we increment s. Figure 3.2 compares the 



theoretically predicted boundary points with the experimentally deter- 
mined ones. Unfortunately, determinants of matrices much larger than 
15 X 15 are not easily calculable on typical home computer hardware. 
Out in the wild, minor expansion tends to get the nod over Gaussian 
elimination when working with symbolic matrices with many variables. 
The analysis above suggests that minor expansion is well suited to han- 
dle many variables, and its advantage only increases when matrices are 
sparse. Multiplications in minor expansion always deal directly with 
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the entries of the matrix. In contrast, the fact that Gaussian ehmina- 
tion involves adding terms to entries means that easy-to-multiply-by 
polynomials (e.g., 0, integers, monomials) not in the first row can be 
"corrupted" after an elimination step. 

It is worth noting that a Gaussian elimination variant introduced by 
Lee and Saunders [5] takes advantage of Os to eliminate unnecessary di- 
vision. However, is not the only easy-to-multiply-by polynomial. We 
hope that this style of cost estimation could lead to better prediction 
of when Gaussian elimination should be used over minor expansion, 
though this prediction depends on implementation details of the algo- 
rithms. 

4. Optimizing minor expansion by row permutation 
Minor expansion is asymmetrical, passing through rows from top 



to bottom in the implementation given in Section 2.1, which we con- 
sider in this section. We mentioned in Section [3] that minor expansion 
takes advantage of easy-to-multiply-by polynomials. Our goal in this 
section is to milk as much advantage out of these easy-to-multiply-by 
polynomials as possible. 

It would be helpful to first precisely define what "easy-to-multiply- 
by" means. Making the same simplifying assumptions about multi- 



plication as in Subsection 3.1 (namely, that calculating pq requires 
T{p)T{q) integer operations, each of which has the same cost), we can 
calculate the cost of minor expansion on A, which we call Cm(^), as 

Cm{A) = J2 5^T(a|,|,)T(det(sub(A [\J\ - 1], J \ {j}))). 

JC[n] jeJ 

Note that, although | det(y4)| is invariant under row swaps and trans- 
position of A, this is not in general the case for Cm{A). However, in 
practice, calculating the cost for all 2n\ row and column permutations 
would take an amount of time comparable to that of the determinant 
calculation itself. 

It is therefore desirable to estimate the cost contribution of each 
entry or row — without taking into account other entries or rows — of 
a matrix and use that to sort the rows, reducing the sorting opera- 
tion to O(n^) time which, in practice, is negligible compared to the 
determinant calculation time. 

The cost contribution of an entry depends on: 

Number of terms: The cost of calculating pq is T{p)T{q). If we 
were only performing a single multiplication, then this would be 
the only important attribute with regards to cost contribution. 
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The cost contribution of a row depends on: 

A function of the number of terms of each of its entries: 

The sum, sum of squares, number of nonzero entries, and so on. 
Number of hnearly independent terms: If, for example, the 
i X i minors all have the same form as polynomials and the 
{i + 1)*^ row has k linearly independent terms, then calculating 
the (i + 1) X (i + 1) minors results in, at most, k terms per term 
of the ? X z minors. 

Recall that absolute value of the determinant of a matrix is invariant 
under row swaps. Given a cost contribution estimate for each row, 
sorting rows by descending cost contribution generally speeds up minor 
expansion. 

To see why this is, consider three polynomials p, q and r. The cost in 
integer multiplications of calculating pqr depends on order of calcula- 
tion. For instance, {pq)r takes T{p)T{q) + T{pq)T{r) integer multipli- 
cations. The more consolidation of like terms occurs when calculating 
pq, the further T{pq) is from its maximum of T{p)T{q). A lower cost 
results from choosing the first two polynomials to have as much consol- 
idation of like terms when multiplied as possible and choosing the third 
to have as many terms as possible. (The second goal is less obvious 
from these equations. It is clearly useful in the case where there is no 
consolidation of like terms, and it must be balanced with the first goal 
outside of this case.) 

Although the cost contribution of rows is not yet well-defined (as, 
indeed, they cannot be without considering all other rows in a matrix), 
experiment shows that even these vague ideas are useful in practice. 



Our findings, shown in Figure 4A, are that sorting strategies are most 
effective when about half of a matrix's entries are 0, where sorting cut 
off, on average, approximately | of the calculation time. We suspect 
that this is a reflection of the increased variance in cost contribution 
of rows. 
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